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Abstract 


Dynamic responses of PEM fuel cells are crucial for mobile applications such as in automobiles. There are four main transient processes in a 
PEM fuel cell, namely, species transport, electric double layer charge/discharge, membrane hydration/dehydration, and heat transfer. In this study, 
a rigorous transient model has been developed, accounting for all four transient mechanisms. The dynamic characteristics have been analyzed, 
corresponding to various changes in working conditions, such as relative humidity and/or cell voltage. Moreover, by using three different membrane 
types, Nafion 112, Nafion 115, and Nafion 117, the effect of membrane thickness on the cell dynamic performance has been investigated, and the 
importance of heat transfer effects on the cell dynamic responses have been highlighted. 


© 2006 Elsevier B.V. All rights reserved. 
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1. Introduction 


With the increasing concerns of sustainable energy and 
environmental issues, polymer electrolyte membrane fuel cells 
(PEMFCs) are attracting more attention from the public and gov- 
ernments. The implementation of PEMFCs, replacing traditional 
combustion engines in automobiles, promises a potential clean- 
energy future. In the past decade, significant improvements have 
been achieved and trial buses, driven by PEMFCs, are in suc- 
cessful operation around the world. However, many technical 
barriers still remain to be resolved before fuel cell power sys- 
tems can be extensively commercialized. For example, the freeze 
startup of PEMFCs at subzero temperature conditions is a com- 
plex problem for automobile applications. It entails three-phase 
modeling (in the gas diffusion layer), namely ice, liquid water 
and water vapour. Only a complete study involving all three 
phases is able to predict timescales of freeze startups as well as 
the negative impacts on reliability and durability of subzero oper- 
ation. In this study, only a single phase is considered, providing 
the building blocks for future multi-phase flow modeling. 

In our previous work [1], a two-dimensional isothermal 
model is developed to investigate the transient phenomena in a 
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fully humidified PEMFC. In the fully humidified case, the mem- 
brane hydration process is of no concern, and the cell dynamic 
performance is determined by the transport of the reactant gas 
species and the electrical double layer charge/discharge process 
only. Our analysis indicates that the time scale for the electrical 
double layer charge/discharge process is on the order of 10 ws; 
the time scale of species transport is longer, but still within 2 s, 
even at severely flooded conditions. Therefore, the cell presents 
excellent dynamic characteristics in the fully humidified case. 
However, experimental observation in our laboratory suggests 
that the dynamic response time in an actual operating cell is 
much longer, and is of the order of minutes, perhaps due to the 
combined effect of slow membrane water uptake, liquid water 
formation in the GDL, and heat transfer processes. Hence, a 
more comprehensive transient model is required to account for 
the water transport and heat transfer processes. 

The water transport in PEM fuel cells entails two compet- 
ing effects: on one side, the membrane electrical conductivity 
increases with the membrane water content; thus the membrane 
should be as hydrated as possible to facilitate the proton trans- 
port. On the other side, however, excess water may accumulate 
in the cathode gas pore and catalyst layer, and block the oxygen 
transport, which causes a significant concentration overpotential 
loss. A dynamic water balance is hence necessary for a sound 
operation of the cell. In general, two modeling approaches of 
membrane water transport exist: hydraulic models and diffusive 
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Nomenclature 


a 


water activity; density of active area in the catalyst 
layer (m? m~?) 

molar concentration (mol m7?) 

specific heat (Jkg~! K) 

hydraulic diameter (m) 

mass diffusivity of species (m? s7!) 
Faraday’s constant 96487 (C mol`!) 

mass transfer coefficient (m s7!) 

heat transfer coefficient (W m~? K) 
reference exchange current density (A m7?) 
phase current density (A m~?) 

thermal conductivity (W m`! K) 

number of electrons transferred in the half cell 
reaction; unit normal vector 
electro-osmotic drag coefficient (H20/H*) 
Nusselt number 

pressure (atm) 

heat source (J m~?) 

universal gas constant 8.314 (J mol! K) 
reaction rate (A m7?) 

relative humidity 

source terms, entropy (J K7!) 

Sherwood’s number 

time (s) 

temperature (K) 


Greek letters 


ee 


transfer coefficient 

porosity 

overpotential (V) 

membrane water content 

density (kg m~?) 

electrical conductivity (S m7!) 

equivalent molecular weight of dry membrane 
(1.1 kg mol7!) 

relative humidity 


Subscripts and superscripts 


anode 

activation 

cathode 

effective value 

gas phase 

the ith species 

membrane; mass transfer 
source term of species 
ohmic 

reference state 

reversible 

solid phase; species transport 
saturation pressure (atm) 
source term of temperature 
water vapor 

membrane water 


b source term of charge 
0 standard conditions (273 K and 1 atm); inlet con- 
ditions 


models. In the hydraulic model, the convective flow of liquid 
water caused by a pressure gradient is included. However, exper- 
imental work [2] has shown that the application of a pressure 
difference between the cathode and anode did not have a large 
effect upon the drag coefficient, which is contrary to what is 
suggested by the hydraulic model. Hence, the diffusive model 
is more frequently employed in the literature. Springer et al. [3] 
were one of the pioneers working with the diffusive membrane 
hydration mechanisms. They implemented a one-dimensional 
isothermal model which accounts for the variation of membrane 
water content. Empirical correlations between the water activity 
and membrane water content are developed based on experimen- 
tal data measured at 30°C. Other transport properties such as 
the water diffusivity and membrane electrical conductivity are 
considered as functions of membrane water content. Indepen- 
dently, Yang [4] developed an empirical expression that related 
the membrane water content with water activity based on their 
experimental results, which correspond to the results of Springer 
et al. Kulikovsky [5] used a new set of correlations based on the 
experimental work of Hinatsu et al. [6] which were measured at 
80°C, and implemented a model similar to Springer et al., which 
differentiates the water transport in the membrane from the water 
vapor transport in the GDL and gas channel. In the catalyst layer 
(CL), two phases are coupled together through an equilibrium 
assumption and the PEM water concentration is converted to 
the water vapor concentration through mathematical relations. 
Using a similar mathematical technique, Um and Wang [7] pro- 
posed a single-domain water transport model which converts the 
water concentration in the electrolyte to water vapor concentra- 
tion both in the catalyst layer and the membrane layer, and the 
resulting water transport equation can be casted into a general 
form which is valid for the whole domain. This approach is also 
employed in the present study and a more detailed description 
of this approach will be given in the next section. It should be 
noted that all of the above mentioned models assume that the 
hydronium (or proton) concentration in the membrane is con- 
stant, and only the water concentration varies spatially in the 
membrane. Berg et al. [8] established a model which accounted 
for the hydronium concentration variation in the membrane. Sev- 
eral fitted parameters are explored to simplify their model. Most 
recently, Baschuk and Li [9] derived a more general membrane 
model. They considered the membrane as being composed of a 
mixture of solid matrix, hydronium ions and liquid water, and 
described the transport phenomena in the membrane using the 
Stephen—Maxwell equation. 

A heat transfer sub-model is also necessary for a rigorous 
transient analysis. Compared to water transport, heat transfer in 
the cell is much easier to model and, hence, there is less discrep- 
ancy in the literature. However, almost all physical parameters 
and transport processes are related to temperature. To reduce the 
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computational complexities of the system, heat transfer effects 
are neglected by most of the researchers. Ju et al. [10] pre- 
sented a single phase non-isothermal model which coupled the 
heat transfer with electrochemical reactions and mass transport. 
A parametric study on the GDL thermal conductivity is con- 
ducted and the heat release in each part of the cell is analyzed 
in detail. The authors conclude that the thermal effects become 
more critical at higher current density and/or lower GDL ther- 
mal conductivity. In the work of Rowe and Li [11], the heat 
transfer in the cell is considered as the combination of the fol- 
lowing two mechanisms: conductive heat transfer in the solid 
matrix; and convective heat transfer in the pores, with local 
thermodynamic equilibrium between the two phases. The heat 
generation/adsorption effects during the water phase change are 
also taken into account in this work. Berning et al. [12] later pre- 
sented a similar conductive-convective heat transfer model. In 
their model, the temperature in solid and fluid phases are solved 
independently, and a relatively arbitrary heat transfer coefficient 
is used to account for the heat exchange between the two phases. 

Most of the above-mentioned studies treated the water and 
thermal management separately and none of them dealt with 
the transient phenomena of PEMFCs. Reviewing the literature, 
we found that most of the studies are focused on steady state 
modeling and very few of them concentrated on transient phe- 
nomena. Berg et al. [13] presented a transient discharge model 
for the cathode, and used it for parameter tuning and deter- 
mination of liquid water in the membrane electrode assembly 
(MEA). Stumper et al. [14] used a similar approach to deter- 
mine the MEA resistance and electrode diffusion coefficient. 
Natarajan and Nguyen [15] presented a two-dimensional two- 
phase transient model for the cathode. Both multi-species flow 
and capillary flow of liquid water are accounted for in their 
model. Wang and Wang [16] developed a three-dimensional, 
single-phase model for the whole cell; the evolution of water 
accumulation and the corresponding dynamic responses of the 
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cell performance are studied. However, all above-mentioned 
transient models are isothermal and none of them investigated 
the cell as a water transport and thermally coupled system. 

Therefore, the objective of the present study is to develop a 
comprehensive transient model which incorporates water trans- 
port and heat transfer in the cell, and to investigate the effect of 
dynamically changing operating condition, such as a dynamic 
change in cell potential or relative humidity, on the cell perfor- 
mance. 


2. Model formulation 


In this study, a two-dimensional modeling domain, identical 
to our previous study [1], is adopted and shown in Fig. 1. The 
convective flow in the gas flow channel is not considered in the 
present model since the flow direction is normal to the 2D mod- 
eling plane. The detailed description of mass and charge transfer 
processes were illustrated in [1]. Therefore, only those matters 
related to water transport and heat transfer will be elucidated in 
this study. 


2.1. Water transport 


In the present study, a pseudo two-phase water transport 
model similar to Springer and Kulikovsky’s approach is used. 
That is, adsorbed water is present in the electrolyte and vapor in 
the GDL. In the catalyst layer, both phases coexist in local ther- 
modynamical equilibrium. The mechanisms of water transport 
are schematically shown in Fig. 1(a). 

In the GDL, the phase change of water is neglected and water 
exists as vapor only. Furthermore, due to the two-dimensional 
nature of the current model (cross-plane model), convection 
between the gas channel and the catalyst layer is also ignored, 
with the convection at the interface between gas flow channel and 
GDL taken care of via mass transfer coefficients in the bound- 
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Fig. 1. (a) Schematic diagram of water transport and heat transfer in a PEM fuel cell; (b) modeling domain and boundary conditions. 
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ary conditions (see further below). Hence, the transient process 
of water transport in the backing layer is governed by diffusion 
only: 


dc. 
sa = V(DEVew) (1) 


where cw is the water vapor concentration, £g is the porosity of 
the GDL, and D$, is the gaseous water diffusivity. 

In the membrane, the movement of protons tends to drag 
water molecules with them from the anode towards the cathode, 
the so-called electro-osmotic drag effect. On the other hand, 
water is produced at the cathode catalyst layer. The electro- 
osmotic drag is thus confronted by back diffusion of produced 
water, which points from cathode to anode. The governing equa- 
tion for water transport in the membrane can be stated as: 


dCwm 


ot 


na> 
em = V(DEV cwm) — V (Jn) 2) 
where Cwm is the water concentration in the membrane, £m is the 
volume fraction of electrolyte in the membrane, Df) is the liq- 
uid water diffusivity in the membrane, nq is the electro-osmotic 
drag coefficient, and Jm is the protonic current density in the 
electrolyte phase. 

In the catalyst layer, both ionomer water and water vapor 
coexist and are taken to be in thermodynamic equilibrium 
locally. We assume water diffusion in the electrolyte and water 
vapor diffusion in the pores of the catalyst layer. The governing 
equation for the equilibrium system is written as follows: 


dc. dc 
es + Em = V(D8Vcw) + V(DYV cwm) 
l R-V (Si ) (3) 
2F ' Po 


where the two terms on the left hand side represent the water 
sorption in the voids and electrolyte, respectively. The third and 
fourth terms on the right hand side represent the water production 
and electro-osmotic drag, respectively. Water is produced at the 
cathode side catalyst layer only, hence the third term on the 
right hand side vanishes at the anode catalyst layer. The water 
concentration in the membrane, Cwm, is defined as: 


PmA 

Cwm = Sy (4) 
where pm is the dry membrane density, W is the equivalent 
molecular weight (1.1 kg mol~!) of the dry membrane, and À 
is the membrane water content, which denotes the number of 
water molecules per sulfonic acid group within the membrane. 
As mentioned previously, several correlations for the membrane 
water content have been developed in the literature. Among 
them, the formula developed by Springer et al. [3] was used 
most extensively. However, the experimental data presented by 
Springer et al. was measured at 30°C, which has non-negligible 
differences from the data at 80°C. This can be easily seen in 
Fig. 2. Therefore, an empirical formula for the experimental 
data at 80°C is employed in this study [5]: 


O Experimental data at 30 °C 
o Experimental data at 80 °C 
—— §Springer's correlation 
Kulikovsky's correlation 


Water content, A 


0.00 0.25 0.50 0.75 1.00 
Water activity, a 


Fig. 2. Water uptake of Nafion membrane at equilibrium with water vapor. 


à = 0.3 + 6a[1 — tanh(a — 0.5)] 


43.9Ja i + tanh (a (5) 


In the above equation, the membrane water content is a unique 
function of the water activity, a, which is defined as: 


CwRT 


a= — 
Psat 


(6) 


The water saturation pressure, P*', is determined through 


[3]: 
log) P** = —2.1794 + 0.02953(T — 273) — 9.1837 
x 10-5(T — 273)? + 1.445 x 1077 (T — 273} 
(7) 


where PS" is in units of atm and T, the temperature, is in units 
of Kelvin (K). 

Strictly speaking, Eq. (5) is only valid for a< 1. However, 
it was extrapolated to the region 1 <a<3 in the present study 
to account for the supersaturation conditions which commonly 
appear in single-phase models, such as in Kulikovsky’s work. 

From Eqs. (4)-(6), it can be seen that the water concentra- 
tion in the membrane, Cwm, can be converted to the water vapor 
concentration, cy, through the mathematical relation [5]: 


dewm — pmRT da 
dew | WP» da 


(8) 


Consequently, using the single-domain water transport 
approach [7], Eqs. (1)—(3) can be combined into a general form: 


eff OCw 
Ew 
ot 


1 Nd = 
= V(De'Vew) — Ri -V ( s Ja) (9) 
where the effective porosity (e) and effective diffusivity (D) 


of water transport can be evaluated according to Eq. (8) in each 
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layer. For example, in the catalyst layer it has the form: 


d 
eci = ég + €m ( cm) (10) 
w 
d 
pet = pe + p™ ( =) (11) 
WwW 


Once the membrane water content is obtained from Eq. (5), 
the liquid water diffusivity and electro-osmotic drag coefficient 
can be determined with the empirical correlations given in [5]: 


Be A — 2.5 
D™ = 4.1 x 10 (=) h + tanh (+=) (12) 


1, à <9, 
na = (13) 
0.117. — 0.053, A>9. 


The membrane electrical conductivity, om (S m7 !), is also a 
function of water content and is expressed as [3]: 


1 o4 
= (0.0051392 — 0.00326 12688 | —— -- 
om = ( a | (= T 


(14) 


2.2. Heat transfer 


Similar to the water transport, convection effects on the heat 
transfer are also neglected in our two-dimensional model. The 
general form of the energy equation can then be stated as: 


oT ; 
(pcp) = VRMVT) + Ò (15) 
where (Pcp) represents the effective heat capacity of each 
layer. For example, in the porous GDL we have: 


(pcp) = eC pcp), + (1 — €N(pcp), (16) 


That is, the effective value is determined by the volume frac- 
tion of each phase in the layer. The specific heat (cp) of a Nafion 
membrane is unavailable. The specific heat of a PTFE based 
structure is used instead. Similarly, the effective thermal con- 


catalyst layer it reads: 
kof = ewkw + €mkm + Esks + €gkg (17) 


where éw, €m, €s and €g, represent the volume fraction of liq- 
uid water, membrane electrolyte, solid matrix and the gas pore, 
respectively. Hence, the summation of the € terms equals 1. 

The last term Q in Eq. (15) represents the heat source. The 
total heat generation ina PEM fuel cell is comprised of reversible 
and irreversible heat releases. The irreversible heat release can 
be further divided into activation heat generation and ohmic 
heating, since in the present formulation the activation polariza- 
tion actually includes the overpotential due to the mass transfer. 
Therefore, the total heat generation is: 


Ò = Ove + Ü act + Oonm 


Ri JZ J 
= |—_| (TAS) + |n* Ri| + = + = 18 
= (TAS) + |n Ril aet t gef (18) 
Ů/ act N. 
rev ohm 


where AS represents the entropy change of the overall reac- 
tion, 7 is the activation overpotential, R; is the reaction rate, os 
and om are the electronic and protonic conductivities, respec- 
tively, J, and Jm are the electronic and protonic current densities, 
respectively. Rj, Js and Jm will be defined in the next section. 

It should be noted that the first two terms — the reversible and 
activation heat generation — are related to electrochemical reac- 
tions and, hence, only valid in the catalyst layers. The protonic 
current flows in the electrolyte. Therefore, the protonic ohmic 
heating source terms are considered only in the membrane and 
the two catalyst layers. Similarly, the electronic ohmic heating 
term is considered in the bipolar plates, electrodes and catalyst 
layers, with the membrane excluded. 


2.3. Governing equations 


The species and charge transport equations have been already 
formulated in our previous study [1]. The governing equations 
for the whole system can now be summarized as follows (the 
source terms in the subsequent equations are summarized in 


ductivity, kf, is evaluated for each layer. For example, in the Table 1): 
Table 1 
Source terms for the conservation equations in each region 
Sk So ST 
Bipolar plate 0 0 a 
Electrode 0 0 J 
as 
1 J2 J2 
Anode CL Sm = ap Ri Ri In Ril + oop + ger 
i s m 
Smo = —V (44m) 
1 Ri , L pn 
Cathode CL So, = — zF Ri Ri =| (TAS) + |n * Ril 4 Set t ef 
nay 1 = 
Smo = -V (Im) — ap Ri 
2. 
Membrane 0 0 i, 
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Ic; ; 
Species : a = V(DËVc;) + Sk (19) 
Electrons: 0 = V(otVV,) — So (20) 
Protons: 0 = Viattiv Vin) + So (21) 
-OT 
Energy : (pcp) = V(k"VT) + Sr (22) 


Here, Eq. (19) describes the diffusion of O2, N2, and H2O 
for the cathode side, and diffusion of H2, and H20 for the anode 
side. Eqs. (19)-(22) are closely coupled through the nonlinear 
Butler—Volmer equation: 


j ! ci N” Qan F (Ve — Va) 
=a ex —— = 
i JO, ref tae p RT s m 


en F 
ap | — V, v|} (23) 


and the current densities are related to the phase potential 
through the following expressions: 


J, = -otf V V, (24) 
Ím = 0V Vm (25) 


2.4. Boundary and initial conditions 


The boundary conditions for water transport and heat transfer 
processes are illustrated in Fig. 1(b). Water transport modeling 
usually involves the implementation of internal boundary con- 
ditions at the membrane/CL interface and/or the liquid/vapor 
interface. However, by using the one-equation approach in this 
study (Eq. (9)), the complexity of applying internal boundary 
conditions is simply bypassed. As a result, a convective mass 
transfer boundary condition is applied at the channel/GDL inter- 
face which needs the water vapor concentration to be specified 
in the gas channel. We use 


= ff 
n(DE V cw) electrode side 
= hy (Cw, 0|channel side — Cwlelectrode side) (26) 


where n denotes the unit normal vector, cy, is the water vapor 
concentration in the gas channel which is calculated based on 
the relative humidity, operating temperature and pressure: 
iS; ps% 

Cw 0o = —— 27 

w,0 RT ( ) 
Here, 3 denotes the relative humidity (RH) which is one of the 
input parameters, P**' is the saturation pressure and determined 
by Eq. (7). The parameter Am in Eq. (26) is a mass transfer coef- 
ficient. The mass transfer of water vapor from the gas channel to 
the electrode surface is analogous to heat transfer in a channel 
with heat flux on one boundary and the other three boundaries 
being insulated. Therefore, hm can be calculated from the Sher- 
wood number correlation for uniform surface mass flux [17]: 


hmdp 


w 


Sh = = 2.712 (28) 


where dp represents the hydraulic diameter of the gas chan- 

nel, treating the gas flow in the channel as laminar. For the 

other boundaries in the water transport region, no-flux conditions 

apply 

OCw 

on 

In the heat transfer domain, convective heat transfer boundary 

conditions are applied on all gas channel boundaries: 


ACV T)| solid side = A'T(To|channel side — T|sotid side) (30) 


where the temperature at the gas channel To is kept constant at 
353 K; hry is the convective heat transfer coefficient, which can 
be derived from the Nusselt number correlation [17] 


hrd 
Nu = = 


=0 (29) 


= 3.61 (31) 
I 

Here, k; denotes the thermal conductivity of the gas mixture, i.e. 

humidified hydrogen at the anode gas channel and humidified 

air at the cathode gas channel. 

For the other boundaries in the heat transfer region, periodic 
(or thermal insulation) boundary conditions apply. To leading 
order and for the purpose of this study, this is the case for fuel cell 
stacks. The initial conditions are simply the steady state fields 
of the previous operating point. 


2.5. Numerical procedures 


The conservation equations are solved using COMSOL 
Multiphysics™, a commercial software package based on finite 
element methods. The GMRES linear system solver with the 
Geometric Multigrid preconditioner is selected to accelerate the 
computation. Variant time stepping is used in which a small con- 
stant time step is used at the first few time steps to capture the 
relatively fast species transport phenomena, and a larger time 
step is chosen for the rest of the time period for the slower water 
transport and heat transfer processes. To accurately describe 
the electrochemical phenomena, the element size in the catalyst 
layer is refined to no more than 1/10 of the catalyst layer thick- 
ness and the resulting mesh has approximately 4 x 10° degrees 
of freedom (DOF). It takes about 1 h to run a transient case on a 
64-bit Linux machine with dual core processor at 2.8 GHz and 
4GB memory. 


3. Results and discussion 


In our base case, the gas feeds are assumed to be fully humid- 
ified (100% relative humidity) at a temperature of 353 K and 
a pressure of 2 atm. Water condensation and the ensuing two- 
phase phenomena exceed the scope of the present study but 
will be incorporated in our future studies. Thus, the relative 
humidity may exceed | in some regions where the water vapor 
can be assumed to be in a supersaturated state. The thermal 
and physical properties used in the model are listed in Table 2; 
other parameters such as the geometric dimensions and kinetic 
parameters can be found in [1]. For comparison purposes, three 
kinds of membranes (Nafion 112, Nafion 115, Nafion 117) are 
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Table 2 
Thermal and physical properties of some materials @ 353 K [18,22] 


Density (kg m~’) 


Thermal conductivity (W m~! K) Specific heat (Jkg~! K) 


Bipolar plate 1900 
Electrode support (PTFE) 2200 
Catalyst layer 2100 
Membrane (Nafion) 1980 
Hydrogen 0.069 
Air 0.9950 
Water vapor 0.632 
Liquid water 972 


21 710 
1.3 (Ballard AvCarb® P150) 1050 
0.8725 1050 
0.445 1050 (PTFE) 
0.204 14400 
0.03 1010 
0.023 1960 
0.67 4197 


modeled in this work. Studies [19] have revealed that Nafion 
115 and Nafion 117 membranes have very similar microstruc- 
tures, but the Nafion 112 microstructure differs from that of 
Nafion 115 and Nafion 117 in the crystalline component of the 
PTFE backbone. Consequently, Nafion 112 possesses different 
membrane properties from Nafion 115 and Nafion 117, such as 
the permeability, water uptake values, etc. In this study, how- 
ever, we assume all Nafion membranes have the same physical 
properties. They differ from each other only by the thickness. 
The thickness of a dry membrane Nafion 112, Nafion 115, and 
Nafion 117 is 50.8 um, 127 um, and 177.8 um, respectively. 
They increase linearly with an increase in membrane water 
content à [3]. Therefore, the membrane thickness in PEMFCs 
under the working conditions is usually adjusted by using certain 
expansion coefficients to account for the swelling effect during 
the membrane water uptake. In this study, some typical thick- 
ness values which are widely employed in the literature (50 um 
for Nafion 112 [20], 164 um for Nafion 115 [21], and 230 pm 
for Nafion 117 [12]) are used to qualitatively demonstrate the 
effect of membrane thickness on the cell dynamic performance. 

Proton transport in the membrane will cause membrane 
potential losses due to ohmic effects. It is usually called mem- 
brane ohmic overpotential. Physically, the ohmic potential loss 
is proportional to the resistor length and inversely proportional 
to the electrical conductivity. Provided that the conductivity 
of Nafion 112, Nafion 115, and Nafion 117 are the same, it 
can be inferred that the membrane ohmic overpotential will be 
determined by the membrane thickness, and hence the overpo- 
tential increases are in the order of Nafion 112, Nafion 115, and 
Nafion 117. In other words, the cell performance decreases in 
the same order. This conclusion is verified by the polarization 
curves shown in Fig. 3. Moreover, for the same kind of mem- 
brane, such as Nafion 117, the ohmic overpotential is determined 
by the electrical conductivity which, in turn, is determined by 
the membrane water content. The membrane water content, and 
thus the electrical conductivity, is higher under isothermal con- 
ditions where we find generally lower temperature than for the 
non-isothermal case. Therefore, the isothermal case shows bet- 
ter performance than the non-isothermal one, as shown in Fig. 3. 
This trend also applies to Nafion 112 and Nafion 115, although 
they are not shown in the figure. 

To better understand the transient processes within a work- 
ing PEM fuel cell, separate water transport and heat transfer 
studies are conducted using Nafion 117. Fig. 4 exhibits the 
dynamic responses of current density with respect to a step volt- 


age change from 0.5 to 0.6 V. In Fig. 4(a), the water transport 
sub-model is switched off and the membrane electrical conduc- 
tivity is assumed to be constant at 0.10 S cm~! [22]. The figure 
shows that an abrupt undershoot exists at the initial stage which 
is caused by the species transient transport, then the cell gradu- 
ally reaches steady state at around 12 s. Recall that the time scale 
of species transfer is about 1 s, as shown in our previous study 
[1]. Therefore, the prolonged time periods must be brought by 
the heat transfer process. In Fig. 4(b), the water transport sub- 
model is switched on, whereas the cell is assumed isothermal at 
353 K. It can be seen that due to the sluggishness of water trans- 
port, the cell current density varies more smoothly compared to 
the first case, and the state steady is reached at around 40 s. If the 
water transport and heat transfer sub-model are applied simul- 
taneously, as shown in Fig. 4(c), it is seen that the cell current 
density experiences an overshoot right after the undershoot and 
the steady state is further extended to around 80s. The under- 
shoot is still caused by the species transport where the species 
concentration remains low as in the higher current density state. 
As for the overshoot, there are two reasons. Firstly, the tempera- 
ture varies relatively slowly as the working conditions are being 
changed. Therefore, the temperature remains high at the initial 
stage as in the higher current density state. Higher temperature 
results in better performance of the electrochemical reactions, 
thus an overshoot follows. A more important reason for the cur- 
rent overshoot is that the decrease of temperature induces the 
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Fig. 3. Polarization curves at steady state. 
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Fig. 4. Transient variation of current density with respect to a step voltage change (0.5-0.6 V) using membrane Nafion 117: (a) water transport sub-model switched 
off; (b) heat transfer sub-model switched off; (c) both water transport and heat transfer sub-model switched on. 


decrease of saturation pressure as well (refer to Eq. (7)), which 
in turn results in the increase of water activity and membrane 
water content (refer to Eqs. (6) and (5)). The results indicate 
that the membrane water content is slightly increased when the 
cell voltage is switched from 0.5 V to 0.6 V, although there is 
more water being produced at 0.5 V. This again demonstrates 
the significance of the inclusion of thermal effects. 

Fig. 5 compares the dynamic responses of the cell corre- 
sponding to a step voltage decrease from 0.7 to 0.6 V, by using 
three different membrane schemes. For Nafion 117 (Fig. 5(c)), 
the trend of variation is much like the reverse of Fig. 4(c). They 
do have slight differences though. It is observed that the volt- 
age decrease case (Fig. 5(c)) reaches steady state at about 70s, 
slightly faster than the voltage increase case (Fig. 4(c)). The 
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dynamic response of the cell using Nafion 115 is similar to that 
of Nafion 117, but the response time is reduced to around 40s 
due to the reduced membrane thickness, as shown in Fig. 5(b). 
For Nafion 112, water can be readily removed or supplied in 
the thin membrane and Fig. 5(a) shows that there is no obvious 
undershoot before the cell reaches steady state at about 15s. 
The membrane water content at two membrane cross sections 
(line 1 & 2, shown in Fig. 1(b)) is shown in Fig. 6. Line | and 2 
represent the region of the membrane under the flow channel and 
the channel land, respectively. At constant cell voltage of 0.6 V, it 
is seen that the water content of Nafion 117 is higher than Nafion 
112 in both locations. This is due to the temperature and electro- 
osmotic drag effect. From the polarization curves in Fig. 3, it 
follows that the current density for Nafion 112 and Nafion 117 


20 
Time (s) 


40 60 20 40 60 80 


Time (s) 


(c) 


Fig. 5. Transient variation of current density with respect to a step voltage change (0.7—0.6 V) using membrane: (a) Nafion 112; (b) Nafion 115; (c) Nafion 117. 


240 H. Wu et al. / Journal of Power Sources 165 (2007) 232-243 


12 T 7 + - 
N —4— Line1, Nafion 117 
re —s— Line2, Nafion 117 | 
=A- Linet, Nafion 112 
=O- Line2, Nafion 112 
10 
e 
= 
g 9 
= 
fe) 
O 
S 8 
© 
z 
Fd 
6 
5 1 1 A f 
0 0.2 0.4 0.6 0.8 1 
Cathode Anode 
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at 0.6 V are about 0.95 and 0.54 Acm7?, respectively. There- 
fore, the temperature and electro-osmotic drag in the Nafion 
112 membrane is higher than in Nafion 117. The water content 
and electrical conductivity in Nafion 112 are hence lower than 
in Nafion 117. Furthermore, it is apparent that the water content 
always decreases from the cathode side towards the anode. In 
Nafion 117, the two curves intersect at around 0.75 of the non- 
dimensionlized membrane thickness, which suggests that at this 
location the membrane water content under the channel area is 
the same as under the land. Past this point, the water content 
under the channel becomes higher. In contrast, in a Nafion 112 
membrane the water content under the land is always higher than 
that under the channel. 

Fig. 7 demonstrates the transient variation of the water con- 
tent along the center cross section of Nafion 117 (line 3) and 
along the membrane/catalyst layer interface (line 4) (shown in 
Fig. 1(b)) with respect to a relative humidity decrease from 1 
to 0.5. Fig. 7(a) shows that the membrane water content varies 
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Fig. 8. Transient variation of the membrane water content corresponding to a 
relative humidity increase from 0.5 to 1 (Nafion 117, cell voltage 0.6 V): (a) on 
line 4; (b) on line 3 (refer to Fig. 1(b)). 
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faster under the channel than under the land, because the water 
vapor under the channel area is more readily removed. The varia- 
tion diminishes at around 50 s under the channel, but diminishes 
at around 120s under the land area. In the membrane transverse 
direction (line 3), Fig. 7(b) exhibits that the water content first 
decreases at the two interfaces (anode and cathode), then the 
variation evolves gradually towards the center. At around 2s, 
the water content in the whole region begins to drop. The water 
content distribution follows a parabolic-like curve during the 
variation. The evolution subsides at around 120s, similar to the 
dynamics shown on line 4. 

Fig. 8 presents the transient variation of the water content 
on line 3 and line 4 (refer to Fig. 1(b)) of Nafion 117 with 
respect to a relative humidity increase from 0.5 to 1. On line 
4 (Fig. 8(a)), the water content under the flow channel adjusts 
to the relative humidity increase very quickly due to the fast 
water vapor diffusion from the gas channel. At around | s, the 
water content is almost uniform. Between 1 and 10s, the water 
content in the whole membrane increases rather uniformly. At 
around 40s, the membrane region under the channel reaches 
steady state first; the water content under the land region keeps 
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Fig. 9. Transient variation of the temperature corresponding to a relative humid- 
ity decrease from 1 to 0.5 (Nafion 117, cell voltage 0.6 V): (a) on line 4; (b) on 
line 3 (refer to Fig. 1(b)). 


on increasing and reaches steady state at about 60s. Fig. 8(b) 
represents the transient variation of the water content on line 3. 
Very similar to Fig. 7(b), the variation starts from the two ends, 
then evolves towards the middle. At around 2s, the water con- 
tent in the whole region begins to increase in a parabolic manner. 
The evolution slows down at 60s and the cell reaches equilib- 
rium. Comparing Fig. 8 with Fig. 7, we can also find that the cell 
responses faster when the relative humidity increases rather than 
decreases, last but not least because of nonlinear water diffusion 
in the membrane. 

The temperature variation in Nafion 117, corresponding to a 
relative humidity decrease (1—0.5), is presented in Fig. 9. Gener- 
ally speaking, the variations of temperature along both directions 
(line 3 and line 4) are more uniform compared with the variation 
of the water content, and the transient variation curves are almost 
parallel to each other. The thermal conductivity of hydrogen is 
several times higher than that of air (refer to Table 2). There- 
fore, more heat is removed from the anode gas channel and the 
temperature decreases from the cathode to the anode, as shown 
in Fig. 9(b). Also, the reversible cathodic reaction is exothermic 
while the reversible anodic reaction is endothermic, i.e. more 
heat is released in the cathodic reaction than in the anodic reac- 
tion. The steady state of temperature distribution is achieved 
around 120s, the same as for the water content. 

Fig. 10 compares the dynamic responses of three different 
membranes for the relative humidity changes. It can be seen 
that the response time increases from Nafion 112 to Nafion 117, 
due to the increased membrane thickness. Moreover, it is seen 
that the cell experiences longer response time when the relative 
humidity decreases rather than increases. Using Nafion 117 as an 
example, the response time is about 120 s as the relative humidity 
is decreased from 1 to 0.5. However, the response time is only 
about 60s as the relative humidity is increased from 0.5 to 1. 
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Fig. 10. Transient variation of the current density with respect to the relative 
humidity changes (solid line: 1-0.5; dashed line: 0.5-1): (a) Nafion 112; (b) 
Nafion 115; (c) Nafion 117. 
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The cell dynamic performance under low relative humidity 
working conditions is also studied and the results for Nafion 117 
are demonstrated in Fig. 11. Itclearly shows that the cell dynamic 
characteristics at lower relative humidity are totally different 
from that at fully humidified operation (Figs. 4(c) and 5(c)). 
This is mainly because the membrane can take up more water 
under lower relative humidity conditions than at fully humidi- 
fied conditions. In Fig. 11(a), the cell voltage is switched from 
0.5 to 0.6 V. An undershoot occurs at the initial stage due to 
the low species concentrations. With the recovery of reactant 
species supply and due to the relatively higher temperature and 
hence better electrochemical performance at the initial stage, the 
cell current density increases and an overshoot appears. On the 
other hand, a large amount of water is accumulated in the mem- 
brane at the initial higher current density state. Therefore, the 
membrane dehydrates gradually due to the reduced water pro- 
duction at lower current density. Accordingly, the cell current 
density decreases gradually till it reaches an equilibrium state at 
around 160s. In Fig. 11(b), a reverse transient curve exhibits a 
cell voltage decrease from 0.6 to 0.5 V. An undershoot follows 
an overshoot due to the species transient, then the current den- 
sity gradually increases during the membrane hydration process 
and reaches steady state around 140s. Comparing Fig. 11 with 
Figs. 4(c) and 5(c), it surprisingly reveals that the response times 
are almost doubled when the relative humidity is reduced by a 
half. 


4. Conclusions 


A comprehensive transient model which accounts for all four 
transient processes, namely, species transport, electric double 
layer charge/discharge, membrane water uptake, and heat trans- 
fer, is implemented in this study. By using different membrane 


types, i.e. Nafion 112, Nafion 115, and Nafion 117, the effect 
of the membrane thickness on the dynamic cell performance 
is studied. The results show that the cell dynamic response 
time increases with the membrane thickness and the Nafion 112 
membrane presents the best dynamic performance. 

The membrane hydration/dehydration processes are analyzed 
through the transient water content distribution along two dif- 
ferent directions. The results reveal that the membrane becomes 
hydrated/dehydrated faster in the region under the channel than 
under the land. The inclusion of heat transfer processes has sig- 
nificant influence on the dynamic cell response. The current 
density varies monotonically if the cell is isothermal. In contrast, 
both current overshoot and undershoot occur if the temperature 
effect is accounted for. 

Numerical studies on the change of cell working conditions 
show that the dynamic cell response is not symmetric under 
symmetric changes in working conditions. The cell response is 
slower if the cell voltage is increased or the relative humidity 
is decreased, both resulting in a current density decrease. On 
the other hand, when the cell operates at low relative humidity 
conditions, the dryer membrane can uptake and reserve more 
water than the fully humidified membrane. Thus, the dynamic 
response of the cell becomes more sluggish. The response time 
can be doubled if the relative humidity is reduced by half. 

In summary, heat transfer effects have a substantial impact 
on transient modeling results. This might be amplified in our 
future two-phase modeling efforts. 
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